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ABSTRACT 

We have developed a numerical simulation code that treats the transport and accel- 
eration of charged particles crossing an idealized oblique, non-relativistic shock within 
the framework of pitch angle transport using a finite-difference method. We consider 
two applications: 1) to study the steady-state acceleration of energetic particles at an 
oblique shock, and 2) to explain observed precursors of Forbush decreases of galactic 
cosmic rays before the arrival of an interplanetary shock induced by solar activity. For 
the former, we find that there is a jump in the particle intensity at the shock, which is 
stronger for more oblique shocks. Detailed pitch angle distributions are also presented. 
The simple model of a Forbush decrease explains the key features of observed precur- 
sors, an enhanced diurnal anisotropy extending several mean free paths upstream of the 
shock and a depletion of particles in a narrow loss cone at ~ 0.1 mean free path from 
the shock. Such precursors have practical applications for space weather prediction. 

Subject headings: acceleration of particles - - shock waves - - cosmic rays - - solar- 
terrestrial relations 
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1. Introduction 

An amazing variety of astro-physical phenomena 
can be attributed to magnetohydrodynamic shocks, 
and in particular to their ability to accelerate charge d 
particles to high energies ( Blandford fc Eichler 1987 ). 



Shock acceleration is believed to account for the bulk 



of the galacti c cosmic rays (Axford 1981 ) and solar 
cosmic rays ( Lee fc Ryan 1986 ) observed near the 
Earth. There is a variety of mechanisms by which 
shocks can accelerate particles, as indicated by di- 
rect observations within the solar system. Diffusive 



shock acceleration (Krymskii 1977, Axford, Leer, & 
skadron 1977. Bell 1978), based on the first-order 



Fermi acceleration mechanism (Fermi 1954, Parker 



1958), can accurately account for acceleration at the 



Earth's bow shock (Eichler 1981, Ellison, Mobius, & 
Paschmann 1990) and traveling interplanetary shocks 
(Lee 1983, Kennel et al. 1986] ), whereas stochastic ac- 
celeration, or second-order Fermi acceleration (Fermi 
1949), is apparently the dominant mechanism in the 



vicinity of cometary bow shocks ( Fcrasawa & Scholcr 
1989|). 



Several lines of evidence indicate the distributed 
acceleration of energetic ions over a wide range of he- 
liolongitudes during gradual solar flare/coronal mass 



ejection (CME) events ( 


Mason, Gloeckler, & Hoves- 


tadt 1984, Lee & Ryan 198£, 


Reames 1990 
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it is not clear whether the first- or second-order Fermi 
mechanism is responsible for initial particle accelera- 
tion near the Sun. As CME shocks move outward 
from the Sun, they become traveling interplanetary 
shocks, and as such a shock passes a detector, in- 
tensity enhancements or anisotropy variations of en- 
ergetic particles are often recorded. Evidence that 
such intensity enhancements can represent true accel- 
eration, not mere pile-up, comes from major spectral 
changes (|McKibben 1972|, [Reames et al. 1997[). At en- 



ergies < 1 MeV, Lhun 



l Liulus 



lu is Lluai uviut'iice LliaL pail 
can 1 be dLLeleialed uuL uf Lhu sulai wind pupulaliuii 



(Gosling et al. 1981), presumably by diffusive shock 
acceleration. However, for energies > 10 MeV, all 
observations to date of ionic charge states near the 
times of interplanetary shock arrivals are inconsistent 
with a solar wind origin, indicating that the ions were 



originally accelerated out of coronal material (Boberg 



Tylka, fe Adams 1996| ). In addition, Tan et al. (1989) 
provided compositional evidence that ions of ~ 1 MeV 
observed at the time of shock passage represent the 



same population as ambient solar energetic particles. 
The above evidence indicates that particles are in- 
jected into the interplanetary medium while the shock 
is still near the Sun, and later this particle popula- 
tion may be further affected by the shock as it propa- 
gates through the inner heliosphere. This leads one to 
the question of how an oblique shock (i.e., a shock in 
which the magnetic field is neither parallel nor per- 
pendicular to the shock normal) further accelerates 
existing populations of energetic charged particles. 

Note that close to an oblique shock, the diffusion 
approximation does not provide an accurate descrip- 
tion of the spatial or directional distribution of en- 
ergetic particles. The goal of the present work is to 
examine the spatial and pitch angle distribution that 
arises due to the transport and acceleration of an ex- 
isting particle population near an oblique shock. To 
the author's knowledge, this represents the first solu- 
tion of the pitch angle transport equation on both 
sides of an oblique, non-relativistic shock for non- 
ultrarelativistic particles. (For ultrarelativistic parti- 
cles at oblique shocks, for which one can set E ss pc, 
see Kirk fc Heavens 1989] ; for parallel shocks an d 
non-relativistic particles, see Kirk fc Schneider 1989.) 



Therefore, we have started with the simplest case 
of an oblique shock with a constant magnetic field 
on either side in a medium with a spatially uniform 
scattering mean free path (Figure 1). It is hoped 
that this will lay the groundwork for further stud- 
ies of particle acceleration integrated into the frame- 
work of a pitch angle transport equation that will 
consider other effects for particular types of shocks; 
these could include a realistic magnetic field configu- 
ration, a spatially dependent scattering amplitude, a 
self-consistent treatment of wave generation and pitch 
angle scattering, or other effects. 

The basic process of charged particle acceleration 
by planar, parallel shocks has been worked out in the 
diffusion approximati on (e.g., Krymskii 1977| , [Bland- 
ford fc Ostrikcr 1978| ). Decker & Vlahos (1986) and 



Jokipii (1987) found that for oblique shocks the rate 
of particle acceleration increases with the shock an- 
gle, i.e., the angle between the magnetic field and the 
shock normal. In addition, there is a characteristic 
length over which the particle distribution increases 
upstream of the shock, given by D/u — vA||/(3it), 
where D is the coefficient of spatial diffusion (due 
to pitch angle scattering), u is the fluid speed rela- 
tive to the shock, v is the particle speed, and A|| is 
the spatial mean free path along the magnetic field. 
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the rate of particle acceleration vs. the magnetic 



AU (Palmer 1982|), and for u ~ 0.002c and v = 0.1c 


field-shock angle ( 


Decker & Vlahos 1986 




Lieu et al. 


(E = 94 MeV), this length scale is on the order of 


1994 




Naito & Takahara 1995 


), the injection of parti- 



1 TA || , or 1.4 to 5 AU. Thus a model with a constant 
diffusion coefficient does not explain the more local- 



cles from a thermal population and its decreased effi- 



ciency for larger magnetic field-shock angles (Baring 



ized increases observed as an interplanetary shock 



passes the Earth. One explanation is that the dif- 
fusion coefficient may decrease near a shock due to 
increased magnetic turbulence generated by accclcr- 



Ellison, & Jones 1993), compositional selection effects 



(Ellison, Jones, & Eichler 1981), and the relative im- 
portance of first- and second-order Fermi acceleration 



ated particles (e.g., McKenzie fc Westphal 1968 , Lcc 
198^, iReamcs, Barbier, fc Ng 1996|) . On the other 
hand, the large changes in particle anisotropics that 
are sometimes observed during shock passage (e.g. , 



Evenson, Meyer, fc Yanagita 1982| , [Heras et al. 1995] ) 



indicate a breakdown in the diffusion approximation, 
so a detailed treatment should also consider the pitch 
andle distribution. Recently, the EPAM instrument 
on the ACE spacecraft has been able to measure com- 
plcx pitch angle distributions close to the time of 
shejek passage (Armstrong & Holland 1998); mod- 
els of shock acceleration should attempt to reproduce 
such distributions. Furthermore, changes in the pitch 



at shocks (Krulls & Achterberg 1994). 

Here we adopt a different approach by incorpo- 
rating the treatment of an oblique shock disconti- 
nuity into a finite-difference numerical solution of a 
pitch angle transport equation. Previously, numeri- 
cal techniques have been developed to treat the pitch 
angle dependent trans port of energetic particles in 
the inner helio s phcrc (|Ng fc Wong 1979j |Schlutci| 
1985| , |Earl 19871 , [Ruffolo 199l| , |Pauls fc Burger 1994|) , 



recently including t he effects of adiabatic decelera- 
tion and convection (Ruffolo 1995 , Hatzky 1996 , Lario| 
1997 , Kocharov et al. 199S ), and these have been 



anj ic distribution of rclativistic ions (mainly gaiac- 
tic cosmic ray protons) as an interplanetary shock 
api roaches, representing precursors to a b 'or bush de- 



used to analyze observations of solar energetic par- 
ticles (Bieber et al. 1980 , Bicber, Evenson, fc Pomcr- 



antz 1986, Droge, Wibberenz, fc Klecker 1990, Kallen- 



crease accompanying the shock ( Forbush 1938 ~ Berry 



rode, Wibberenz, & Hucke 1992, Ruffolo, Khum- 
lumlcrt, fc Youngdcc 1998|) and solar neutron decay 



particles (Ruffolo 1991,~ |Droge, Ruffolo, fc Klecker 
1996| ). Including the injection of particles from other 



& Hess I942~b'orbush & Lange 1942 ) can be measured 
by ground- based neutron monitors and could provide 



advance warning of approaching shocks (Bieber & 
Evenson 1997), which can induce geomagnetic storms 



and affect satellites, communications, and power grids 
at Earth. Clearly the implementation of such a warn- 
ing system would rely on an accurate understanding 
of the pitch angle distribution ot energetic ions up- 



sources, such as a traveling interplanetary shock, 
has proven more difficult. Previous authors have 
set the simulation boundary at the shock and as- 
sumed an ad hoc particle injection function, say, 
Q, finding the function that best fits the observed 



data for a given event and energy range (Heras et 
al. 1995L iLario 1997L iKallcnrodc & Wibberenz 19971 



stream ot an interplanetary shock. 

Kirk & Schneider (1987a) presented an approx- 
imate analytic solution of a pitch angle transport 
equation for highly rclativistic particles near a par- 
allel shock discontinuity (i.e., one in which the mag- 
netic field is parallel to the shock normal), giving 
insight into the variety of length scales correspond- 
ing to higher orders of anisotropy near the shock 
(the length scales also apply to non-relativistic par- 
ticles). This work was extended to oblique shocks 
by Kirk & Heavens (1989). Monte Carlo techniques 
have provided important information on the pitch 



Schnci- 



dcr 1987b 




Naito & Takahara 1995|), the final spec- 


trum (Kirk & Schneider 19871, Ballard & Heavens 


19Snl Lieu 


ct al. 1994 


, Naito & Takahara 1995 


)• 



Kallenrode 1997a, b). In principle this allows one to 
examine the dependence of Q on physical parameters 
of the shock, though there are a number of such inter- 
related parameters, which complicates the interpreta- 
tion of the association between Q and any individual 
parameter. Our approach is in a sense the reverse: we 
examine what parameter dependence should be ex- 
pected based on certain physical processes. Since dif- 
fusive shock acceleration is a manifestation of particle 
transport in the vicinity of the shock, we include both 
sides of the shock discontinuity in the calculation, 
with special treatment of particles crossing the shock 
based on the conservation of the magnetic moment, 
to examine both the transport and further accelera- 
tion of energetic charged particles near the shock. We 
are able to examine the effects of the shock-field an- 
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gle and the form of the pitch angle scattering on the 



de Hoffmann- Teller (shock) frame (dc Hoffmann & 



steady-state particle distribution in space and pitch 



Teller 1950) where the electric field is zero. The en- 



angle for an idealized model of an oblique shock. As 
a further application of our method, we apply it to 
explain the observed precursors of Forbush decreases 
before the arrival of an interplanetary shock. 

2. Methodology 
2.1. Overview 

Before describing the methodology in detail, we 

prpapnt. a hripf nvprvipw nf sl-inpk arpplprat.inn TVip 



tire energy change due to mechanisms 1 and 2 is ac- 
counted for by transforming the particle momentum 
from the local fluid frame into the shock frame, con- 
sidering the energy-conserving shock encounter in the 
shock frame, and then transforming the momentum 
into the new local fluid frame (Decker 1983). [Note, 
however, that there is still a lateral drift along the 
shock, which is important, e.g., for the acceleration 
of anomalous cosmic rays at the solar wind termi- 



nation shock (Pesses, Eichler, fc Jokipii 1981, Cum- 



prnlrpss nf "Hiffiisivp slinrV arpplpratinn" rpfprs t.n tlio 



mings, Stone, fc Webber 1985 ).] For an oblique shock 



acceleration of charged particles as they are repeat- 
edly scattered back and forth across a shock front, 
which for the present purposes is a discontinuity in 
the fluid speed and magnetic field strength. This pro- 
cess does not require the diffusion approximation and 
can be examined by lower-level descriptions in terms 
of pitch angle scattering ( Kirk fc Schneider 1987a) or 
particle orbits in a disordered magnetic field (Decker 
fc Vlahos 19851) . 

Originally the literature discussed two mechanisms 
of acceleration at shocks: 

1. In the first-order Fermi acceleration mechanism 



(Fc fini 1954 -P-a ricor 1958 ] ), particloo gain energy when 



sca ttering off of converging magnetic field irrogulari 

ties. At a shock, this occurs because the field irregu- 
larities are convected with the fluid flow (in the regime 
in which the Alfven speed is smaller than the fluid 
speed relative to the shock), and the upstream fluid 
flows toward the shock faster than the downstream 
flui d flows awa} r . — In this regime, scattering tends to 



isoii ropigo the momentum while preserving its magni 
tude in the local fluid frame. When a particle heads 
toward the shock (in the local fluid frame) in either 
direction, the frame transformation leads to a higher 
momentum in the new local fluid frame, which is con- 
served until the particle encounters the shock again. 
After two such crossings of the shock, the particle has 
a higher momentum in its original reference frame. 
2. In the shock drift mechanism (Schatzman 1963), 



particles drift along an oblique shock front due to the 
sharp gradient in the magnetic field, and this drift is 
along the direction of the electric field so that parti- 
cles can gain a substantial amount of energy in one 
encounter with the shock. 

More recently, it has been shown that the distinc- 
tion between these two mechanisms vanishes in the 



the term diffusive shock acceleration is now generally 
used to refer to this unified description of particle ac- 
celeration that includes both mechanisms 1 and 2. 

The basic mechanism of particle acceleration at 
an oblique shock, as described above, includes two 
processes: a) frame transformations, as in first-order 
Fermi acceleration, and b) a change in pitch angle 
as the particle encounters the shock (which does not 
occur for a parallel shock). Processes a) and b) are 
related to parallel and perpendicular changes in the 
fluid velocity, respectively. Thus they are directly 
analogous to the two components of adiabatic de- 
celeration in the solar wind (e.g 



Webb & Gleeson 



1979, Ruffolo 1995), which are associated with the di- 



vergence of the wind parallel and perpendicular to 
the magnetic field and have been termed "an in- 
verse Fermi effect" and "betatron deceleration," re- 
spectively. Just as betatron deceleration is an effect 
of adiabatic focusing (also known as magnetic mir- 
roring) when viewed in the local fluid frame (Ruffolc 



1995| ) , which in turn represents the conservation of the 



magnetic moment p ± /(2meB), for the case of shock 
acceleration the change in pitch angle as particles en- 
counter the shock approximately conserves the mag- 
netic moment as well (see §2.5). If we assume that 
the magnetic moment is actually conserved, then we 
do not need to concern ourselves with the details of 
the particle orbit near the shock; we treat the en- 
tire particle-shock interaction as a single event, and 
consider whether a particle ultimately crosses or is 
reflected by the shock, which is assumed to depend 
only on the initial pitch angle. 

Figure 2 provides a schematic illustration (for a 
non-relativistic particle) of acceleration at an oblique 
shock. [For an analogous illustration for adiabatic de- 
celeration, see Figure 1 of Ruffolo (1995).] The com- 
ponents of the particle velocity, vu and v±, refer to 
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motion parallel and perpendicular to the magnetic 
field, respectively (the pitch angle is the angle with 
respect to the positive vn axis). Let us take the up- 
stream direction to be to the right; then the upstream 
and downstream flow speeds with respect to the shock 
frame are negative ("U" and "D," respectively, in Fig- 
ure 2). A particle with < is moving from the up- 
stream side toward the downstream side. The veloc- 
ity relative to point "D" is greater than that relative 
to point "U," so if the particle crosses to the down- 
stream side, the magnitude of its velocity is greater 
in the new local fluid frame, and this velocity magni- 
tude is conserved by further scattering in that frame. 
Similarly, a particle with vu > can cross from down- 
stream to upstream, and again the magnitude of the 
velocity increases in the new local fluid frame. A third 
type of shock encounter (for oblique shocks) is when 
a particle from upstream (v\\ < 0) reflects back up- 
stream with i>|| — > —vu in the shock frame; again, we 
see that the magnitude of the velocity increases in the 
fluid frame. 

In addition to the frame transformations described 
above, i.e., standard first-order Fermi acceleration for 
parallel shocks, at oblique shocks there is the addi- 
tional process of a change in pitch angle as a particle 
crosses the shock. At first glance it is not immedi- 
ately obvious why this process should lead to acceler- 
ation as opposed to deceleration. To see this, consider 
again a particle with v\\ < crossing the shock from 
upstream to downstream (Figure 2). In the shock 
frame the large-scale magnetic field configuration is 
static, so the particle speed in this frame, + vj_, 

is conserved. If we assume approximate conservation 
of the magnetic moment, then v± will increase, since 
the magnetic field is stronger on the downstream side. 
As can be seen in Figure 2, the speed relative to the 
downstream fluid ("D") increases; the dotted line in- 
dicates a curve of constant speed in the downstream 
fluid frame. Thus this process leads to an additional 
increase in energy in the new local fluid frame, be- 
yond that due to the frame transformation. Similarly, 
a particle with > is moving from downstream to 
upstream, and approximate conservation of the mag- 
netic moment implies a decrease in vj_. Then the 
speed relative to the upstream fluid ("U") increases; 
the dashed line indicates a curve of constant speed in 
the upstream fluid frame. The change in pitch angle 
can be quite large, so this can lead to a substantial in- 
crease in the energy gain at oblique shocks. There can 
also be a significant effect on the pitch angle distribu- 



tion, which makes particle transport and acceleration 
at oblique shocks more difficult to model than that at 
parallel shocks. 

2.2. Transport Equation 

To our knowledge, this work represents the first 
treatment of oblique shock acceleration within a finite 
difference simulation of pitch angle transport, so as a 
first step this work considers only a simple, planar, 
oblique shock with straight magnetic field lines on 
either side (Figure 1). We consider only subluminal 
shocks so that we can work in the de Hoffmann- Teller 
reference frame (which we also refer to as the "shock 
frame" ) in which the shock is stationary, the fluid flow 
is paral lel to the magnetic field, and the electric field 
is zero (|de Hoffmann fe Teller 195C| ). 

Ruffolo (1995) provided an equation of focused 
pitch angle transport (eq. [11] of that paper) that 
included solar wind effects, such as adiabatic decel- 
eration and convection, to first order in u/c in an 
Archimedean spiral magnetic field. When we consider 
a region of constant magnetic field on either side of 
the shock, that equation reduces to 



dF(t,fj.,z,p) 
dt 



-—fivF{t,fi,z,p) 
az 
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d <p(n) d / uv 
dfi 2 d^i V ' c 2 



F(t,(i,z,p), 



(1) 



where F = d 3 N/(dzdfidp) is the density of 

particles in a given magnetic flux tube, 
t is the time in the shock frame, 
/i is the pitch angle cosine in the wind 
frame, 

z is the distance from the shock along the 
magnetic field in the shock frame, 

p is the particle momentum in the wind 
frame, 

v is the particle velocity in the wind frame, 
u is the solar wind speed along the 

magnetic field in the shock frame, and 
ip is the pitch angle scattering coefficient. 

In equation (|l]), the first term on the right hand 
side represents the effect of streaming, the second 
is for convection (including the relativistic correction 
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for transforming the streaming speed into the shock 
frame), and the third is for pitch angle scattering. 
The particle density in a flux tube, F, is related to 
the phase space density, /, by F = 2irap 2 f, where 
a(z) oc \/B(z) is the cross-sectional area of the flux 
tube. We use F in the simulations, following Ng and 
Wong (1979), because we can easily design the nu- 
merical finite difference method to strictly conserve 
this quantity during streaming and convection, even 
when a(z) and B(z) are spatially varying. The pitch 
angle scattering coefficient, <p([i), is expressed as 



^) = A\nr\\-^) 



(2) 



where A is the scattering amplitude and q controls 
the form of the scattering coefficient. This expression, 
originally derived in the context of quasi-linear theory 
(Jokipii 1971), is adopted as a convenient and widely 
used parameterization. In this work we will consider 
5 = 1, for isotropic scattering, and q = 1.5, which is in 
the range of 1.3 to 1.7 inferred by Bieber et al. (1986) 
for interplanetary scattering. Roughly speaking, q > 
1 implies a deficit in scattering near a pitch angle 
of 90° (fj, = 0); the effect of such a deficit on shock 
acceleration has also been examined by Kirk (1988). 

2.3. Eigenfunction Expansion 

In addition to our numerical solution of equa- 
tion (|l|) (see §2.6), it is possible to find analytic solu- 
tions in the steady state if we restrict the z domain 
to one side of the shock. Setting dF/dt = 0, we have 

dF A 

dz 



2\pv(l — /iuv/c 2 ) + u] 



9 , 

V 1 



uv\ 
M^j F. (3) 



Note that F is defined in a mixed frame, i.e., for \x 
and p defined in the solar wind frame and z and t 
in the shock frame. To first order in u/c, the analo- 
gous quantity defined in terms of variables in the solar 
wind frame is given by F w = (1 — /j,uv/c 2 )F (Webb 
Sz Gleeson 1979), so we can simplify equation (|3|) to 
read 



dF w 
dz 



1 a -(l- M 2 ) | Mr i^k ) (4) 



2 [/it> + it] dfi dfi 

again neglecting terms of order (u/c) 2 . For a given 
p, this can be solved by separation of variables, 
F w (fi,z) = M(fi)Z(z), yielding Z oc e kz and 



d ,„ 9m m i dM ( u\ 



where a = 2kv/A is an eigenvalue of the equation. 
To avoid divergence as z — > ±oo, we must have 
k < (> 0) for z > (< 0). Following Kirk & Schnei- 
der (1987a), who treated the case of q = 1, M (fi) can 
be expanded in terms of normalized Legendre polyno- 
mials. For q and u/v values of interest, we truncated 
the expansion to n terms and evaluated eigenvalues 
and eigenfunctions using the Mathematica software 
package (Wolfram Research, Inc.). For q = 1, n = 12 
was sufficient to obtain eigenvalues with a relative ac- 
curacy of 10 -5 . For q = 1.5, obtaining an accuracy 
of about 2% for the first positive eigenvalue required 
n rs 80; for other eigenvalues fewer terms were re- 
quired. 

The eigenvalues a closest to zero (corresponding to 
long spatial scales) are shown in Figure 3 for q = 1 
and 1.5 and a range of u/v values. For u/v — 0, the 
numerical values of a are 0, ±14.53, ±42.05, ±83.30, 
±138.3, ...for q = 1 and 0, ±11.3, ±33.0, ±65.6, 
±109, . . . for q = 1.5. For a given q the relative mag- 
nitudes of these eigenvalues are similar, but not iden- 
tical, to those of the Legendre polynomials, £(£+1) for 
£ = 0,1,2,.... The eigenvalues shown in Figure 3 for 
q = 1.5 are about 20% lower in magnitude than those 
for q = 1, with the exception that a\ is lower by a fac- 
tor of about 2.4. Aside from a.\ , the eigenvalues shown 
here are approximately given by on = a°(l ± (3u/v), 
where a? is the value for u/v — 0, "±" follows the 
sign of and (3 is a constant between 2.2 and 2.4. 
Sample normalized eigenfunctions for these eigenval- 
ues are shown in Figure |for q = 1 and u/v = -0.075, 
with particles traveling away from the shock for fi < 
downstream and u > upstream; each eigenfunction 
M (fi) is plotted with the sign for which (M) > 0. In 
cases considered in this paper, uv/c 2 -C 1, so F w w F 
and M essentially gives the angular dependence of F. 

For u/v > 0, there is a zero eigenvalue, eta, with a 
corresponding constant eigenfunction, and a first pos- 
itive eigenvalue, ai, which tends to zero as u/v ^ 0. 
For q = 1, ai ~ 6u/v (with slight deviations for 
higher u/v), implying that k w 3uA/v 2 . Note that in 
the diffusion approximation, the spatial diffusion co- 
efficie nt is given by the classic formula (Jokipii 1966, 
1968, Hassclmann fc Wibbcrcnz 1968 ), 



D 



^l-M 2 ) 1 



dfi, 



(6) 



2 Jo ¥>0«) 

and for the form of </?(//) specified in equation (0), we 



(5) 
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have the well-known expression 
D ^ 1 



A (2-g)(4-g) 



(7) 



For q = 1 we have I? = i> 2 /(3A), and the approx- 
imate z-dependence corresponding to the first pos- 
itive eigenvalue is Z oc e uz / D ; such a dependence is 
also found for the case of q = 1.5. Note that this solu- 
tion only applies to the region where uz < 0, i.e., the 
upstream region. In general, solutions correspond- 
ing to positive (negative) a are valid in the upstream 
(downstream) region, and the constant eigenfunction 
corresponding to «o = is valid in either region. 

The complete solution for F in the steady state 
on either side of the shock is a linear combination 
of the separable solutions. Therefore, the solution far 
upstream of the shock is a superposition of a constant 
and the above solution proportional to e uz / D . Since 

|ck a | 3> eti, the only valid solution far downstream of 

the shock is a constant. Such spatial behavior of F at 
long distances from the shock is in accordance with 
standard results based on the diffusion approximation 
(e.g., Krymskii 1977 ; see also §2.5), though there are 
slight deviations at higher u/v (for q = 1, otx > 6u/v 
by «2% at u/v = 0.1). At shorter distances there 
may be contributions from solutions corresponding to 
other eigenvalues. In terms of the particle mean free 
path, 

3D V 3 (8) 



A 



) A(2-q)(A- q y 
such solutions decay over a distance scale of 



1 

k 



2 v 
a A 



2(2-g)(4-g) 
3a 



A. 



(9) 



(The symbol A refers to the mean free path parallel 
to the magnetic field, A||, unless otherwise specified.) 
For q = 1 and u/v = 0.02, deviations from the solu- 
tion given by the standard diffusive approximation are 
expected upstream within 1/k — 0.13A of the shock 
(corresponding to oiq) and downstream within 0.14A 
(corresponding to a_i); for q — 1.5 those distance 
scales become 0.07A and 0.08A, respectively. 

Note that this analysis alone cannot give the rel- 
ative amplitudes of different solutions; for that, one 
must consider how particles cross the shock. Kirk & 
Schneider (1987a) presented eigenvalues and eigen- 
functions of equation (|^) for the case of q = 1, 
and for the case of ultrarelativistic particle energies 
(E « pc) and a parallel shock, were also able to match 



the upstream and downstream solutions and analyti- 
cally evaluate the steady-state distribution function, 
as well as the power-law index of the resulting spec- 
trum. The analogous calculation for oblique shocks 
was performed by Kirk & Heavens (1989). For par- 
allel shocks, Kirk & Schneider (1989) relaxed the as- 
sumption of ultrarelativistic particle velocities in cal- 
culating the steady-state particle distribution as a 
function of position, pitch angle, and momentum. In 
this work, we treat the crossing of moderately rela- 
tivistic particles across an oblique shock in numerical 
(finite difference) solutions of the pitch angle trans- 
port equation, and we also compare our steady-state 
numerical results with pitch angle distributions and 
spatial decay constants, k, that are expected based on 
the eigenfunctions and eigenvalues of equation ^) , re- 
spectively. 

In summary, the results for steady-state shock ac- 
celeration can be expressed in terms of a superposi- 
tion of separable solutions of equation (3). A numer- 
ical solution of equation (1), in our case by a finite 
difference method, is required in order to determine 
the relative amplitudes of the separable solutions. We 
stress that the numerical method does not explicitly 
involve the eigenvalues and eigenfunctions, so we can 
test the numerical code by verifying that the results 
are consistent with a linear combination of separable 
solutions. 

2.4. Fluid and Magnetic Field Parameters 

Based on equations (133)-(135) and (124) of dc 
Hoffmann & Teller (1950), the jump conditions of a 
non-relativistic shock in the de Hoffmann- Teller frame 



are 



3u| 
5 u n 



2u n 



-u 2 s + -u 2 n (l + t 2 ) 



u A , 



= 



= 



= o, 



(10) 



where a bracketed quantity refers to the difference be- 
tween that quantity on either side of the shock, u s — 
\J (5/3)(p/p) is the speed of sound, UAn = B n j 
is the Alfven speed corresponding to B n (hence the 
Alfven speed is ua — uaz sec 9), and t = tan (9, where 
9 is the angle between the magnetic field and the 
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shock normal. In the above, we used the conser- 
vation of mass flux and assumed an ideal gas with 
an adiabatic index of 5/3. Analogous conditions for 
relativistic shocks were derived by Heavens & Drury 
(1988). 

Equations ( |l0| ) relate u n , u s , ua 7 and t upstream 
and downstream of the shock, so given any four 
of those variables the remaining four can be deter- 
mined. We use the subscript 1 (2) to refer to up- 
stream (downstream) variables. Solutions were ob- 
tained for Usi — uax — 50 km s _1 , typical of near- 
Earth interplanetary space, and selected values of t\ 
and u\ n — U2n- Note that in the high speed limit 
uin/u2n —> 4, i.e., (7 + I) / (~f — 1) for an adiabatic 
index 7 = 5/3. In this limit, the magnetic field ten- 
sion becomes negligible, and the tangential compo- 
nent of the velocity u n t, is nearly conserved, so we 
have £2/^1 - > 4 as well. 

The only fluid parameters that directly affect the 
particle transport are the fluid speed upstream and 
downstream and the magnetic compression ratio, B2I 'B±. 
Since the normal component of B is conserved, we 
have B 2 /Bi = sec 9 2 / sec 9\. Figure || shows the 
shock compression ratio, ux n /u2 n , and the magnetic 
compression ratio, B2/B1, as a function of Au n = 
uin — U2n for selected values of t\ = tan#i. Since 
ti/ti — * 4 in the limit of high Mach numbers, the 
asymptotic value of the magnetic compression ratio 
is + 16i?)/(l + t\), which is less than 4. To ex- 
amine the effect of B2/B1 on cosmic ray transport 
and acceleration near the shock, we perform simula- 
tions for Ux n — U2n = 400 km s -1 and t± = 0, 1, and 
4, corresponding to upstream shock angles of 0°, 45°, 
and 76° and downstream shock angles of 0°, 76°, and 
86°, respectively. Complete fluid and magnetic field 
parameters for these configurations are given in Ta- 
ble 1. We do not consider t > 4 because we expect 
that our model would be inaccurate for nearly per- 
pendicular shocks in that it neglects particle diffusion 
perpendicular to the magnetic field, which is impor- 



tant for highly oblique shocks (e.g., Jones, Jokipii, & 
Ba ling 1998] ). 



2.5. Boundary Conditions 

Next we consider changes in the momentum and 
pitch angle as a charged particle encounters the shock, 
i.e., a sudden change in the magnetic field and fluid 
speed. Following Decker (1983), we assume that in 
the shock frame both the particle momentum and the 



first adiabatic invariant, i.e., the magnetic moment 
p\j '{2meB), are conserved. The above assumption 
is tantamount to considering the adiabatic limit, in 
which the shock is a region where the pitch angle 
changes gradually. The opposite limit would consider 
an infinitely thin shock front, and in this case the 
pitch angle distribution can only be determined by 
computing particle orbits. For highly oblique shocks 
(8 ;> 80°), there is essentially no difference between 
the resulting pitch angle distributions in the two lim- 
its ( Tcrasawa 1979 ). Therefore, the approximation we 
are using should be reliable for highly oblique shocks, 
and at least represents a well-defined limit for less 
oblique shocks. 

The original simulation code of Ruffolo (1995) is 
capable of treating multiple values of the momen- 
tum in order to examine the particle distribution as 
a function of the momentum. However, for simplicity 
we have only treated one value of the momentum in 
this work, and have assumed that the dependence of 
F(t, jj,, z,p) on p is given by a power law, F cx p 1 . 
When particles are accelerated during an encounter 
with the shock, we consider them to be advected from 
a lower value of the momentum to the momentum of 
interest. Therefore, we should evaluate the accuracy 
of the assumption that F has the same dependence 
on fi and z at the lower values of p from which parti- 
cles are advected. For non-relativistic particles with 
u < v, the largest fractional momentum change that 
can be achieved (that for reflection back upstream) is 
< 2ui/v (see Figure 2). In Figure 3 we see that the 
dependence of the eigenvalues a on u/v is linear or 
weaker for u/v values of interest. As will be described 
in §3.1, we do in fact see a systematic error in some 
results that increases with u\/v. For ultrarelativis- 



tic particle speeds (e.g., Kirk & Schneider 1987a) the 
assumption of a power-law dependence in the steady 
state is not a problem because v ss c so that u/v docs 
not vary with momentum. 

Now we turn to the boundary conditions at the 
edges of the simulation region. At a reasonable 
distance from the shock front, we can employ the 
diffusion approximation, in which equation (Q) can 
be expressed as a diffusion-convection equation. In 
the manner of Earl (1974), we set F(t,p, fi, z) ps 
F (t,p, z) +Fi(p,n,z), where F Q = Fx is an 

odd function of fx, and higher-order even functions 
are neglected. By separating the transport equation 
into odd and even parts, integrating each over /x, and 
neglecting the time dependence of Fx (i.e., assuming 
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\dFi/dt\ < \vdF /dz\), we obtain 



dF 
dt 



OS* 
dz 



where 



S z = uF - D 



dF 

dz 



(11) 



(12) 



is the z-flux and the spatial diffusion coefficient, D, 
is given by equation (||). Equation (|T^) should be 
obeyed on either side of the shock, with special jump 
conditions at the shock as described above. Note that 
either F or the phase space density, /, may be used 
in these equations and in the following discussion. 

We note further that in terms of F±(p, /x, z), we 
have 



u 1 



1 



2 2 
V C 



F 



liFidy., (13) 



and if F\ — Sfx, so that <5 is the dipole anisotropy in 
the mixed frame, we obtain 



S z =uF + Z(5-^)f,, 



(14) 



where 5 — uv/c 2 represents the dipole anisotropy of 
F w , the distribution function defined entirely in terms 
of independent variables in the solar wind frame. 

In the diffusion approximation solution of Krym- 
skii (1977), there are boundary conditions that F re- 
main finite as z — > ±oo. In the steady state, equa- 
tion (|l|) requires that S z be constant on either side 
of the shock, so wc have 



Fn = A- 



Sexp(^) 



(15) 



with different constants A and B upstream or down- 
stream. To avoid divergence far downstream, we must 
have B = 0, so Fq is constant, with S z = uFq. Up- 
stream, we have S z = uF u , where F u is Fq far up- 
stream of the shock, representing the existing particle 
population that is accelerated by the shock. 

In the numerical solution of equation (jl]) in terms 
of z we can only treat a finite domain. However, the 
diffusion approximation is valid sufficiently far from 
the shock (see §2.3), with dS z /dz = 0, so setting S z 
at the boundaries of the simulation region is approx- 
imately equivalent to fixing its value at ±oo. At the 
upstream boundary, we set S z to a constant (for a 
given p) that is interpreted as uF u . At the down- 
stream boundary, S z is set to u{F) M at the boundary. 



In the diffusion approximation (eq. [12]), this implies 
that OFq /dz — 0, and also that the anisotropy is given 
by S = uv/c 2 ; the distribution function in the wind 
frame is F w = (1 — [luv/c 2 )F and thus has a dipole 
anisotropy of zero. 

2.6. Numerical Method 

The numerical simulations reported here employed 
the finite difference method of Ruffolo (1995) as sub- 
stantially modified by Nutaro, Riyavong, & Ruffolo 
(in preparation). The latter report will contain de- 
tails of the modifications and testing of the new code. 
Here, for completeness, wc briefly outline the changes 
to the treatment of streaming and convection, and the 
treatment of particles crossing the shock. For a de- 
scription of other aspects of the numerical method 
which remain unchanged, the reader is referred to 
Ruffolo (1995). 

Changes to the treatment of streaming and con- 
vection were motivated by the work of Hatzky (1996), 
who us ed the total v ariation diminishing (TVD) tech- 
nique ( [Bweby 1984 ). While our previous numerical 
method (Ruffolo 1991, 1995) eliminated spatial dif- 
ferencing and hence numerical diffusion in the eval- 
uation of these terms, it required a small step size, 
Az = A^vAt. On the other hand, the TVD algo- 
rithm permits only a very small amount of numerical 
diffusion for significantly larger Az. The result is that 
the new code runs 1 to 2 orders of magnitude faster. 

Another benefit is that convection is now treated in 
each step instead of by occasional jumps as in our pre- 
vious method, yielding much smoother profiles with- 
out spatial averaging. The previous method was ad- 
equate for the transport of solar energetic particles, 
where one could average over a moderate spatial in- 
terval, and the jumps merely yielded small-amplitude 
oscillations in the intensity as a function of time. 
However, for the present simulations of steady-state 
shock acceleration, with multiple reflections from the 
boundaries and a requirement of fine spatial resolu- 
tion, the new treatment of streaming and convection 
was essential in obtaining the smooth profiles shown 
here. The results for steady-state particle acceleration 
at a parallel shock provide a sensitive test confirming 
the accurate treatment of streaming and convection. 

In our implementation, we modified the standard 
TVD algorithm in order to permit a Courant number, 
7 = v z At/ Az, greater than one. This permits an ar- 
bitrary Az, and in particular, for Az = A/ivAt we 
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were able to reproduce results for solar energetic par- 



ticles that used our previous code (Ruffolo & Khum- 
lunEert 1995|). 



The other modification of the code was to treat 
how particles cross the shock, which was also part 
of the streaming/convection step. For each fx-z grid 
point, we determine whether particles will stream/convect 
as far as the shock. If so, we perform a Lorentz 
transformation into the shock frame, allow particles 
to cross or reflect while conserving the magnetic mo- 
ment, and perform another Lorentz transformation 
back into the local wind frame (see §2.1). The TVD 
algorithm effectively splits a cell into fractions of par- 
ticles destined to move to two different spatial loca- 
tions; at the shock we generalize this approach to split 
cells into fractions destined to move to different /i-cells 
as well. The accuracy of this technique is verified by 
the numerical results for steady-state shock acceler- 
ation at oblique shocks, which are consistent with a 
sum of separable solutions of the transport equation 
(§2.3) on either side of the shock. 

3. Results 

3.1. Steady-State Shock Acceleration 

In a steady state, an equilibrium is reached in 
the evolution of the particle distribution function in 
terms of position, pitch angle, and momentum. In 
the present simulations, we assume that the momen- 
tum dependence is given by F oc p~ 7 , so we find 
the value of 7 that yields a steady solution for F in 
terms of z and p. Figure 6 indicates the flux balance 
that determines 7. Well away from the shock, we can 
use the diffusion approximation (see §2.5) to say that 
the net z-flux far downstream is due to convection: 
S z = u%Fq, where F = (F)^. Far upstream, there is 
a balance between convection toward the shock and 
diffusion away, so that S z = U\F U , where F u is the far 
upstream flux. Here we set F u = 0, so there is a net 
outflow of particles from the shock in the downstream 
direction. This is balanced by the p-flux, S p , rep- 
resenting acceleration at the shock of particles from 
lower momenta to the momentum of interest, as well 
as from the momentum of interest to higher momenta, 
for the appropriate steady-state power-law index, 7. 

Before considering oblique shocks, we tested our 
methodology for the case of a parallel shock [Ox = 
62 = 0). Fluid and magnetic field parameters were 
as listed in Table 1. We considered protons with 
speeds v = 0.5 and 0.1c, corresponding to kinetic 



energies of 145 and 4.7 MeV and momenta of 541.7 
and 94.3 MeV c _1 , respectively. We used grid spac- 
ings of Az/X — 0.025 and 0.005, respectively, and 
A/i = 2/95. For convenience, we set vAt = Az. 
Outer boundaries were placed at ±2.5A and ±0.5A for 
v = 0.5c and 0.1c, respectively. The numerical treat- 
ment of the boundary conditions assumed F u = and 
was sufficiently accurate that the proper behavior of 
F in the diffusion approximation was maintained out 
to within a few grid points from the boundaries, and 
the location of the boundaries did not influence F 
near the shock. A full simulation required about 2 
hours of CPU time on a Sun Ultra- 1 workstation. 

Figure 7 (top panel) shows the spatial dependence 
of the pitch angle averaged phase space distribution, 
near a parallel shock for v — 0.5c and q — 1 (z > 

is the upstream region). As explained in §2.2, our 
simulations solve for F ~ d 3 N/(dzdp,dp), the density 
of particles in a flux tube, which is related to / by 
F = 2irap 2 f, where a is the cross-sectional area of a 
flux tube; therefore, / oc BF. For a parallel shock, 
B is equal on the upstream and downstream side, so 
/ oc F . Throughout this section, / is normalized to 

1 far downstream. From Figure 7 we see that for 
a parallel shock, the spatial dependence of is 
simply that expected in the diffusion approximation 
(§2.5), i.e., constant downstream and exponentially 
decaying to zero upstream with k = u/D. The pitch 
angle dependence was also the same as in the diffusion 
approximation. The same results were obtained for 
q = 1.5 and for v = 0.1c with q = 1 and 1.5. 

For a coarser /x-grid spacing, a spurious peak was 
obtained in {f) lll which tends toward zero as A/i — * 
(a remnant of this is barely visible in the top panel of 
Figure 7). Even for a fine grid spacing, one anomaly 
is that a value of 7 = 2.020 was needed for a steady 
state (i.e., to conserve the flux of particles in the sim- 
ulation region) for v/c = 0.5 or 0.1 and q — 1 or 
1.5, whereas acceler ation theory for parallel shocks 
(e.g., Krymskii 1977 ) yields 7 = 3?/2/ (iti — 1*2) + 1 or 
2.034 for these fluid parameters. The assumptions of 
that theory should apply given that /(/i, z) is that ex- 
pected from the diffusion approximation. We believe 
that the systematic error in 7 obtained from the sim- 
ulations arises from the assumption of a power-law 
dependence of F on p, which neglects the momen- 
tum dependence of F(fi, z) (see §2.5). The key prob- 
lem is that the upstream anisotropy of F depends on 
U\/v. As a test, the code was modified to artificially 
add an "extra" anisotropy for lower momenta, i.e., to 
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multiply F upon acceleration by a /z-dependent fac- 
tor to account for the higher anisotropy of F at the 
lower momentum from which a particle was acceler- 
ated, yielding a similar F(fi,z) and 7 = 2.041. We 
conclude that this explanation can in fact account for 
a systematic error in 7 of the observed magnitude, 
and that 7 is more sensitive to the assumption of a 
power-law dependence than is the distribution of par- 



ticl 3S in space or pitch angle. 



Turning to oblique shocks, Figure 7 shows the spa- 
tial dependence of for tan Si = 1 and 4 and for 
q = 1 and 1.5. In all cases, the distribution function 
farther from the shock is consistent with the diffusion 
approximation, with (F)^ constant downstream and 
exponentially decaying upstream with k = u/D. A 
conspicuous feature of Figure 7 is the jump (discon- 
tinuity) in {f) fi at an oblique shock (the finite slope 
is due to the finite grid spacing in z). This feature 
was also found in simulations by Ostrowski (1991), 
Gicseler et al. (1998), and T. Naito (private commu- 
nication, 1998). Gieseler et al. (1998) present a de- 
tailed theoretical and computational analysis of this 
feature, as well as possible observational signatures. 
We find that the jump is stronger for more oblique 
shocks, and weaker for q = 1.5 than for q = 1. The 
amplitude of the jump is on the order of a few percent 
for such fast particles (v = 0.5c), and our simulations 
indicate that the jump is stronger for slower particles 
(v = 0.1c), i.e., a higher u/v. 

Another difference from the case of a parallel shock 
is that for oblique shocks, additional eigenfunctions 
are excited in /(/i, z) near the shock. (If one is not suf- 
ficiently careful in treating the boundary conditions, 
as I was not during the initial stages of this work, 
additional eigenfunctions are also excited near the 
boundaries; discretization errors also yield spurious 
eigenfunctions near the shock, which become negligi- 
ble for 95 /i-grid points as used here.) For all steady- 
state simulations, f(fi, z) was consistent with a sum 
of separable solutions of equation (3). For tan Si = 4 
(61 — 75°), Figure 8 shows the dependence of / on \i 
and z within ±0.8A of the shock, and Figure 9 shows 
/ as a function of /i for z — ±0.05A. (Recall that we 
use n and p to refer to quantities in the local fluid 
frame; thus these plots are for a constant value of the 
local p. A Compton-Getting transformation to the 
shock frame would have no noticeable effect on our 
distribution plots.) For tan Si = 1, the results were 
qualitatively similar but with weaker anisotropies. 

In Figure 9, we see that upstream distributions 



(thick lines) increase with /i up to fi ss 0.7 (with a 
slightly stronger anisotropy than in the far upstream 
region), and for greater fi values, / drops sharply. 
The reason for the sharp drop is that given our as- 
sumption of conservation of the magnetic moment, 
particles with fi > -Jl — B\ / B2 , or 0.85 in this case, 
have come from downstream. A similar drop in / 



has been called a "deficit cone" (Nagashima et al 
1992j ) or "loss cone" effect ( [Bieber k Evenson 1997^ 



for the case of galactic cosmic ray (GCR) depletion 
at high /i upstream of an interplanetary shock, which 
is due to the paucity of GCR coming from down- 
stream (see §3.2). The same effect occurs here be- 
cause the acceleration of particles coming from down- 
stream is weaker than for particles reflected from up- 
stream. The greatest acceleration occurs for parti- 
cles reflected with the greatest change in pitch angle 
(see Figure 2), i.e., for /i slightly below 0.85. Since 
stronger acceleration implies that / is advected from 
lower momenta, and the particle spectrum increases 
with decreasing momentum in this case, the strongest 
acceleration corresponds to the greatest increase in /. 

In the downstream region, particles are redistributed 
in pitch angle because of changes in pitch angle as 
particles cross the shock; the average flux also in- 
creases slightly due to acceleration. It is worth not- 
ing that for a highly oblique shock, most particles 
coming from upstream are in fact reflected, i.e., when 
\fi\ < yl — B1/B2, or in the case of a strong, highly 
oblique shock, for pitch angles more than 30° from the 
magnetic field direction. Another feature of Figures 
8 and 9 is the sharp gradient in / at // = for the 
case of q = 1.5. For this form of the pitch angle diffu- 
sion coefficient, (p(fi) — A\/i\°^(l — fi 2 ) tends to zero 
as fi — * 0. Since the /i-flux, = — (ip/2)(dF/dfi), 
is slowly varying in a near-equilibrium situation, the 
vanishing diffusion coefficient at fi = is able to sus- 
tain an infinite gradient in F at that value. 

We believe that this behavior of / as a function 
of z and fi is not an artifact of the assumption of a 
power-law momentum dependence because when an 
extra anisotropy was artificially added, the fi and z 
dependence (including the jump at z = 0) was not 
significantly affected; this was also the case for paral- 
lel shocks. On the other hand, computed values of 7 
are strongly affected by the power-law assumption, so 
that this code in its present form is essentially unable 
to determine 7. The error in 7 was much weaker for 
v/c = 0.5 than for v/c = 0.1 (because of the lower 
U\/v ratio). As an example, for q — 1 the 7 values 
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required for a steady state with v = 0.5c were 1.965 
and 1.952 for tan6q = 1 and 4, respectively, while 
with v — 0.1c they were 1.985 and 1.787, respectively. 
Otherwise, the results regarding /(/i, z) for v/c= 0.1 
were qualitatively similar to those shown in Figures 
7 to 9 for v/c = 0.5, with much stronger anisotropics 
and jumps in (f)^ at the shock. 

3.2. Precursors of Forbush Decreases 

To demonstrate the versatility of this method, we 
apply it to model Forbush decreases of galactic cos- 
mic rays (GCR) as an interplanetary shock passes the 
Earth (|Forbush 1938] , [Berry fc Hess 1942| , [Forbush fc 
Lange 1942| ), which represent a transient phenomenon 
instead of a steady state. Ground-based neutron mon- 
itors measure secondary neutrons from the impact 
of relativistic, primary charged particles, mainly pro- 
tons, on the upper atmosphere. Due to selective de- 
flection by the Earth's magnetic field, neutron moni- 
tor observations are sensitive to primary cosmic rays 
from specific directions in space, and the worldwide 
network of neutron monitors provides detailed infor- 
mation on their pitch angle distribution, sensitive to 
variations on the order of 0.1%. Precursors to For- 
bush decreases are of practical interest as possible 
predictors of space weather effects on the Earth, such 
as satellite failures, radio fade-outs, power outages, 
etc., several hours or even days before the passage 
of a major interplanetary shock. Several analyses 
of neutron monitor observations have indicated two 
types of precursors to Forbush decreases: 1) an en- 
hanced diurnal anisotropy of GCR, with an excess 
of particles traveling toward the Sun along the inter- 
planetary magnetic field, and 2) a deficit of GCR in 
a "loss cone," i.e., along a narrow range of pitch an- 
gles directed nearly along the interplanetary magnetic 
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We model the Forbush decrease in a rather ide- 
alized manner, assuming the configu ration of Fig- 
ure 1 and neglecting par ticle drifts ( Nishida 1983j . 
Kadokura fc Nishida 1986 ), spatial dependence of the 

spat.t.prinp; mwn frpp path shnrlr r-nrva.t.nrp tno finite 
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focusing, and adiabatic deceleration. Nevertheless, 
we can explain the basic features of the observed pre- 
cursors, verifying their interpretation in terms of par- 
ticle transport in the vicinity of an oblique shock. 



The simulation conditions were inspired by the dra- 
matic CME event of 1997 April 7, which arrived near 
Earth on April 10-11, and for which a possible loss 
cone feature is identified by Bieber & Evenson (1997). 
In this case, the travel time of 3 days indicates a 
shock speed of ~ 600 km s _1 or only 200 km s" 1 
faster than the (typical) solar wind speed. Assuming 
that the shock normal is radial, we take the upstream 
shock-field angle to be the typical "garden-hose" an- 
gle of 45° (tan#i = 1), and as before we assume 
u s i = uai = 50 km s . Thus we find Au n = 133 km 
S _1 , tan#2 = 3.20, and BijB\ — 2.37, which in turn 
implies that particles crossing the shock from down- 
stream have pitch angles aligned with the magnetic 
field to within 40° (p > 0.76). We used q = 1.5, 
which adequately describes interplanetary scattering 
(Bieber et al. 1986). For the upstream boundary con- 
dition, we specify a constant F M (see §2.5), and the ini- 
tial condition sets F to that constant in the upstream 
region and to zero in the downstream region. We 
used A[i = 2/45 (45 //-grid points) and Az/\ = 0.05, 
where A = 0.3 AU. We considered v = 0.75c, corre- 
sponding to a kinetic energy of 480 MeV. For the mo- 
mentum spectrum, we assumed F cx p , according 
to the model proton spectrum of Reinecke, Moraal, & 
McDonald (1996) for the similar polarity solar cycle 
of 1977. However, the simulation results were very 
insensitive to the GCR spectral index. 

Figure 10 shows the omnidirectional GCR inten- 
sity as a function of position (normalized to 1 just 
upstream of the shock). As the shock moves past a 
fixed observer, one sees a gradual precursor decline 
and a slight recovery as the shock approaches. Such 
gradual declines in some observations were noted by 
Cane et al. (1996). In their one figure for such a For- 
bush decrease, that of 1972 Oct 31 (Figure 5 of that 
paper), it is seen that 2 out of 3 neutron monitor sta- 
tions observed a relative peak near the time of shock 
onset. It would be interesting if the omnidirectional 
flux could be estimated from the worldwide neutron 
monitor network for such events for direct comparison 
with the results of numerical simulations. 

At the shock itself, we see a jump reminiscent of 
that found in S3.1 for shock acceleration with F„ = 



(Figu re 7; see also Ostrowski 1991 , Gieseler et al 



1998 ). Thus this model predicts a discontinuous drop 
at the shock followed by a more gradual decline which 
we identify as the declining slope of a Forbush de- 
crease. We note, however, that this model cannot 
hope to accurately reproduce the detailed features of 
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the Forbush decrease itself, given its simplistic as- 
sumptions (see also §4). 

The phase space distribution, /, as a function of fi 
and z is shown in Figure 11, and pitch angle distri- 
butions in the near upstream (z — 0.05A) and far up- 
stream (z = A) regions are shown in Figure 12. We see 
that there is an overall enhanced sunward anisotropy 
both downstream and upstream, where it persists over 
several mean free paths from the shock. Close to the 
shock, we also find the second type of observed pre- 
cursor feature with a sharp decline in the loss cone 
region. 

Our numerical simulations indicate that f(/i, z) up- 
stream is approximately given by a sum of steady- 
state separable solutions. If a steady state has not 
been achieved, separable solutions of the full trans- 
port equation involve an eigenvalue equation similar 
to equation (5) with u/v — > u/v — l/(kvr), where 
r is the decay time. Given the weak sensitivity of 
eigenvalues to u/v (except a\), the longest-lived solu- 
tions should be similar to the steady-state solutions of 
equation (5) with some modification to k. The time 
scale of evolution of a Forbush decrease, ~1 day, is 
much longer than the travel time of particles across 
the distance scales of interest, ~10 minutes, so it is 
reasonable that transient solutions have disappeared 
leaving only nearly stationary solutions of the trans- 
port equation. 

Both types of observed precursors of Forbush de- 
creases are readily understood in terms of the up- 
stream simulation results, as expressed as a super- 
position of separable solutions of equation (3). The 
constant solution in the far upstream region corre- 
sponds to the zero eigenvalue. The contribution of 
the first eigenfunction is opposite in sign to that seen 
in §3.1, corresponding to an anisotropy directed to- 
ward the Sun over a long spatial scale. Turning to 
the next eigenfunction, corresponding to oli, the spa- 
tial decay scale of the steady-state eigenfunction is 
0.074A||, and in the time-dependent simulations this 
feature has a spatial scale length of w 0.08 A || along 
the magnetic field, or ss 0.11A r in the radial direction. 
This eigenfunction appears with the same sign as in 
Figure 4 for shock acceleration (Figure 9) and repre- 
sents a deficit of particles in the loss cone (fj, ;> 0.76) 
that came from the downstream region which is de- 
pleted in GCR, as well as an enhancement of particles 
with [i just below 0.76 which were accelerated during 
reflection from the shock. Several examples of loss- 
cone deficits were given by Nagashima et al. (1994) 



and Sakakibara et al. (1995), and an enhancement for 
certain pitch angles corresponding to reflection and 
acceleration at the shock was reported by Belov et al. 

(1995) . The predicted spatial decay length for these 
features depends only weakly on the shock speed or 
obliquity (i.e., the eigenvalues depend only weakly on 
ui/v, especially for relativistic velocities), but there 
is a significant dependence on q (see eq. [9]). 

The results shown here, with precursor features 
roughly of the magnitude reported by Cane et al. 

(1996) , were for a simulated duration of 7.4 hours, 
which is substantially shorter than the actual time 
it takes a shock to propagate from the Sun to the 
Earth (typically 2 to 4 days). This is probably due 
to our neglect of enhanced spatial diffusion near the 
shock, which must be present in order to account for 
the sharp GCR gradient at the onset of a Forbush 
decrease. Stronger scattering would stem the tide of 
equilibration and lead to a longer duration of these 
features. 

4. Discussion 

The numerical simulation procedure developed in 
this work is the first to solve a pitch angle transport 
equation on both sides of an oblique, non-rclativistic 
shock without assuming an ultrarelativistic particle 
velocity. Here we have shown applications of the 
technique to study oblique shock acceleration in the 
steady state, and to consider the time-dependent 
problem of Forbush decreases and in particular their 
precursors. To explore the capabilities and limita- 
tions of this type of solution, we have started with 
the simplest case of a plane-parallel, oblique shock 
with straight magnetic field lines on either side (Fig- 
ure 1) and a spatially uniform scattering mean free 
path. Clearly this simple model is neglecting a va- 
riety of important processes, which will be discussed 
shortly. 

One limitation in this work was our assumption of a 
power-law dependence of the distribution function on 
the particle momentum. In further work this assump- 
tion should be relaxed, i.e., the simulations should 
treat different values of the momentum, as has been 
done by Kirk & Schneider (1989). This assumption 
has strongly affected calculated values of the parti- 
cle spectrum, but based on test runs that artificially 
compensate for the error, it seems not to have signif- 
icantly affected the spatial and pitch angle distribu- 
tions. This should be checked in future work. Monte 
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Carlo (MC) techniques are undoubtedly more popular 
than the finite difference method used here, and MC 
simulations are very well-suited for exploring certain 
phenomena. Here we use a finite difference method 
for consistency with an existing code that includes a 
variety of effects known to be important for interplan- 
etary transport (Ruffolo 1995) so that the treatment 
of the shock can be readily incorporated into such sim- 
ulations. For transport simulations in general, finite 
difference or eigenfunction expansion techniques have 
some advantages, such as often requiring less com- 
puting time, the absence of statistical error, straight- 
forward extrapolation to reduce discretization error, 
easier treatment of a high dynamic range in particle 
density, and a straightforward and continuous imple- 
mentation of analytic expressions for scattering and 
other proc esses (for a com parison of different tech- 
niques, see Earl et al. 199£ ) . The last of these is also 



a notable limitation in situations where analytic ex- 
pressions are not available, e.g., situations requiring 
the tracing of particle orbits. If the tracing of particle 
orbits is only necessary near the shock, that could be 
implemented in a finite difference code by a "transfer 
matrix" based on a one-time evaluation of the proba- 
bility of various outcomes at the shock for each /i-grid 
point. 



To the author's knowledge, calculated pitch an- 



gle distributions for steady-state acceleration of en- 
ergetic particles at a non-relativistic oblique shock 
have not reported previously. Such pitch angle dis- 
tributions should be directly comparable with recent 
and upcoming in situ observations near interplane- 
tary shocks. As observed distributions are reported 
in greater detail, they can be expected to challenge 
the theoretical models and help indicate which ad- 
ditional physical processes have important effects on 
particle transport and acceleration at interplanetary 
shocks. 

The pitch angle distributions for steady-state shock 
acceleration represent the superposition of pitch angle 
eigenfunctions with different amplitudes, as first de- 
scribed by Kirk & Schneider (1987a). The amplitudes 
of different eigenfunctions can only be determined by 
treating how particles cross the shock. The upstream 
pitch angle distributions obtained here for oblique, 
non-relativistic shocks and mildly relativisitic parti- 
cles (Figure 9) are qualitatively very different from 
those found for parallel, relativistic shocks and ultra- 



case, those authors found that even for a parallel 
shock, eigenfunctions beyond those present in the dif- 
fusion approximation are excited. The upstream pitch 
angle distribution is highly collimated in the direction 
away from the shock. The downstream distribution 
for the parallel, mildly relativistic shock (ui = 0.3c) is 
qualitatively similar to what we find for oblique, non- 
relativistic shocks. Our results can be more directly 
compared with the pitch angle scattering results of 
Naito & Takahara (1995), who treated oblique, mildly 
relativistic shocks with ui = 0.1c-sec6>i. We find qual- 
itatively similar upstream pitch angle distributions 
but very different downstream distributions; for high 
obliquity (and high ui) those of Naito & Takahara 
(1995) are strongly collimated with particles moving 
away from the shock. 

The second application considered here, to pre- 
cursors of Forbush decreases of galactic cosmic rays 
(GCR), demonstrates the ability of this method to 
simulate time-dependent phenomena. The eigenfunc- 
tion analysis alone is sufficient to specify that in a 
steady state, one can have certain types of features 
in the pitch angle distribution over certain distance 
scales. For example, the observed increase in diurnal 
anisotropy well before the onset of some interplane- 
tary shocks (e.g., Cane et al. 1996, Bieber fc Evenson 



1997) is identified with the eigenvalue a\ and a long 



relativistic particles (Kirk & Schneider 1987a, b, Heav- 



ens fc Drury 1988). In contrast to the non-relativistic 



distance scale, and the superimposed excess for \x up 
to « 0.5 and strong deficit thereafter (Nagashima et 
al. 1992, 1994, |Sakakibara et al. 1995] ) is identified 
with the eigenvalue cti and a distance scale of ~ 0.1 A. 
The time-dependent simulations are necessary to ver- 
ify that these modes are in fact excited with the ap- 
propriate sign and a reasonable amplitude, and that 
the time-dependent, upstream distribution is approxi- 
mately represented by a superposition of steady-state 
eigenfunctions (which is not the case downstream). 

We note that simulations of a Forbush decrease 
with a parallel shock did not exhibit features corre- 
sponding to Q!2) which is an example of how not all 
eigenfunctions are excited in every situation. Since a 
localized deficit in a narrow loss cone is in fact not ex- 
pected for a parallel shock, this supports the physical 
explanation that the upstream deficit over a narrow 
range of pitch angles corresponds to particles cross- 
ing fro m the downstream reg ion which is depleted in 
GCR ( [Nagashima et al. 1992| ). 

As stressed earlier, this model of a Forbush de- 
crease is idealized in many ways, though it seems to be 
adequate for describing the key features of upstream 
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precursors, which are of practical interest for short- 
term space weather forecasting (Bieber & Evenson 
by warning of the impending impact of a ma- 
jor interplanetary shock. We are much more hesitant 
to apply this idealized model to the Forbush decrease 
itself. However, according to Wibberenz, Cane, and 
Richardson (1997), the key features of a Forbush de- 
crease can be captured by assuming enhanced scatter- 
ing (a lower mean free path) in the region just down- 
stream of the shock. This could easily be included 
in a simulation technique such as ours and a compar- 
ison with particle distributions observed during the 
course of Forbush decreases could help identify what 
physical processes are crucial to the Forbush decrease 
phenomenon. 

A key motivation for this work is the potential 
to include more physical effects in the future. For 
a more realistic magnetic field configuration, trans- 
port effects such as adiabatic focusing and deceler- 
ation can readily be included, as they have already 
been included in numerical simulations of interplan- 
etary transport. Enhanced scattering near a shock 
almost certainly affects acceleration and transport in 
that region, yet the magnitude and extent of such 
scattering for high-energy particles is not well under- 
stood. Concrete models of the spatial dependence of 
the mean free path near a shock should be developed 
and tested for their success in explaining in situ ob- 
servations. 
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Fig. 1. — Model configuration in the de Hoffmann- 
Teller frame, in which the shock is stationary, the 
fluid flow u is parallel to the magnetic field B, and 
the electric field is zero. 

Fig. 2. — Illustration of changes in pitch angle when 
a particle crosses an oblique shock, a process which 
contributes to particle acceleration. Here uy and v± 
are components of the particle velocity in the shock 
frame. See text for details. 

Fig. 3. — Eigenvalues closest to zero for equation (5) 
as a function of u/v for q = 1 (solid lines) and q — 1.5 
(dashed lines). 

Fig. 4. — Eigenfunctions of equation (5) for eigenval- 
ues shown in Figure 3 for u/v = —0.075 and q = 1. 
For steady-state shock acceleration, the pitch angle 
distribution is a linear combination of such functions. 
Particle motion is directed away from the shock for 
fi < downstream and fi > upstream. 

Fig. 5. — Fluid compression ratio, ui n /u2 n , arid mag- 
netic compression ratio, B2/B1, as a function of the 
normal velocity jump in the shock frame, Au n = 
uin — u 2n , for t\ = tan 6\ = 0, 0.5, 1, 2, and 4. 

Fig. 6. — Schematic illustration of the balance of spa- 
tial and momentum fluxes in diffusive shock acceler- 
ation. 

Fig. 7. — Spatial dependence of the pitch-angle av- 
eraged phase space distribution function, for 
steady-state particle acceleration near a shock (at 
z = 0) for various values of tan 9\ , the tangent of the 
angle between the magnetic field and the shock nor- 
mal, for q = 1 (solid lines) and 1.5 (dotted lines), and 
for v — 0.5c. Farther from the shock, is constant 
downstream and exponentially decays toward zero up- 
stream. The ordinate is normalized to the value far 
downstream. 

Fig. 8. — Phase space distribution of particles as a 
function of \i and z (in units of A) near an oblique 
shock with tan#i = 4 for a) q = 1 and b) q — 1.5. 
Note the changes in the pitch angle distribution near 
the shock (at z = 0). 

Fig. 9. — Phase space distribution of particles as a 
function of /z near an oblique shock with tan 61 = 4 
for q = 1 (solid lines) and q = 1.5 (dashed lines) at 
z = 0.05A (upstream; thick lines) and z = — 0.05A 
(downstream; thin lines). 
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Fig. 10. — Spatial dependence of the pitch-angle aver- 
aged phase space distribution function, (./%, of galac- 
tic cosmic rays near an oblique interplanetary shock 
(tan 0i = 1). As the shock moves (to the right) past 
a fixed observer, one sees a gradual precursor decline 
and a slight recovery, followed by a discontinuous de- 
crease at the onset of a Forbush decrease. 

Fig. 11. — Phase space distribution of particles as a 
function of /i and z (in units of A) near an oblique 
interplanetary shock (tan#i = 1). The sharp drop 
downstream (z < 0) represents the Forbush decrease; 
upstream (z > 0) features represent precursors that 
may be useful for space weather prediction. Superim- 
posed on the enhanced sunward anisotropy, note the 
unusual pitch angle distribution just upstream of the 
shock. 

Fig. 12. — Phase space distribution of galactic cosmic 
rays as a function of fj, near an oblique interplanetary 
shock (tan 0i = 1) at z — 0.05A (solid line) and z = 
A (dashed line), representing near upstream and far 
upstream precursors of a Forbush decrease. 
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Table 1 

Selected Solutions of Shock Jump Conditions' 1 



tan 61 


tan 6 2 




U 2n 


B 2 /B l 








538.0 


138.0 


1.00 


1 


3.88 


541.5 


141.5 


2.83 


4 


15.11 


544.3 


144.3 


3.67 



a For uai — u s i — 50 km s 1 . All velocities 
are in units of km s _1 . 
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u/v 



Sp to higher p 



S^ = u 2 F 



diffusion 



s*=o 



downstream i Convection upstream 



S p from lower p 



1.04 
1.02 
1.00 

0.98 

0.96 
1.04 

1.02 
1.00 

0.98 

0.96 
1.04 

1.02 
1.00 
0.98 
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